Excess mortality

Combining three countries

Data

expected_deaths_monthly <- read_rds("data/expected_deaths_monthly.Rds") %>% 
  mutate(Pandemic = case_when(
    (Year >= 1889 & Year <= 1892) ~ 1890,
    (Year >= 1917 & Year <= 1920) ~ 1918,
    (Year >= 1956 & Year <= 1959) ~ 1957,
    (Year >= 2019 & Year <= 2020) ~ 2020
  ))

Figure 1 - diffs

Monthly

Yearly

Figure 2 - four pandemics

Absolute

focus_years <- c(1890, 1918, 1957, 2020)

# plot_years <- c(seq(1890 - 1, 1890 + 2),
#                 seq(1918 - 1, 1918 + 2),
#                 seq(1957 - 1, 1957 + 2),
#                 seq(2020 - 1, 2020 + 2))

for(country in unique(expected_deaths_monthly$Country)) { 
  
  print(paste0("Plotting ", country))
  
  max <- expected_deaths_monthly %>% 
    filter(Country == country) %>% 
    filter(!is.na(Pandemic)) %>% 
    summarise(Deaths = max(Deaths),
              lpi_flu_hc = max(lpi_flu_hc),
              upi_flu_hc = max(upi_flu_hc)) %>% 
    rowwise() %>% 
    mutate(max = max(c(Deaths, lpi_flu_hc, upi_flu_hc))) %>% 
    ungroup()
  
  print(paste0("Max for x is ", max$max))
  
  for(year in focus_years) {
    
    # skip Spain oldest - no data
    if(year == 1890 & country == "Spain") {
      next
    }
    
    print(paste0("    Year ", year))
    
    plot_data <- expected_deaths_monthly %>% 
      filter(Country == country) %>% 
      filter(Pandemic == year)
    
    name <- paste(country, year, sep = "_")
    
    if (year == 2020) {
      my_title = paste0(country, ": ", 
                        as.character(year - 1), "-", as.character(year))
    } else {
      my_title = paste0(country, ": ", 
                        as.character(year - 1), "-", as.character(year + 2))
    }
  
    # sync breaks across plots
    my_breaks = plot_data$Date[seq(1, nrow(plot_data), 6)]
  
    plot <- plot_data %>% 
      ggplot(aes(x = Date)) +
      geom_rect(xmin = as.Date(paste0(year, "-01-01")),
                xmax = as.Date(paste0(year, "-12-31")),
                ymin = -Inf, ymax = Inf, 
                fill = "#ffffbf", alpha = 0.01, inherit.aes = FALSE) +
      geom_ribbon(aes(ymin = lpi_flu_hc, ymax = upi_flu_hc), alpha = 0.1) + 
      geom_line(aes(y = fit_flu_hc), colour = "grey50") +
      # geom_line(aes(y = median_flu_year), colour = "#2c7bb6") + 
      geom_line(aes(y = Deaths), color = "#d7191c", size = 0.5) + 
      scale_x_date(breaks = my_breaks, date_labels = "%b") +
      xlab("") + 
      # ylab("Observed, 10y median & predicted Deaths") +
      ylab("Observed & predicted Deaths") +
      # labs(caption = "Excludes extreme years") +
      ggtitle(my_title) +
      scale_y_continuous(limits = c(0, max$max),
                         labels = label_number(accuracy = 0.1, suffix = "K", 
                                               scale = 1e-4)) +
      theme_bw() +
      theme(plot.title = element_text(size = 10),
            axis.text.x = element_text(size = 8),
            axis.title.y = element_text(size = 9),
            axis.text.y = element_text(size = 8))
    
    assign(name, plot)
  }
}
[1] "Plotting Switzerland"
[1] "Max for x is 10766"
[1] "    Year 1890"
[1] "    Year 1918"
[1] "    Year 1957"
[1] "    Year 2020"
[1] "Plotting Sweden"
[1] "Max for x is 17278"
[1] "    Year 1890"
[1] "    Year 1918"
[1] "    Year 1957"
[1] "    Year 2020"
[1] "Plotting Spain"
[1] "Max for x is 163422"
[1] "    Year 1918"
[1] "    Year 1957"
[1] "    Year 2020"
plot_spacer() + Spain_1918 + Spain_1957 + Spain_2020 +  
  Sweden_1890 + Sweden_1918 + Sweden_1957 + Sweden_2020 +
  Switzerland_1890 + Switzerland_1918 + Switzerland_1957 + Switzerland_2020 +
  plot_layout(widths = rep(c(2, 2, 2, 1), 3)) +
  plot_layout(ncol = 4) 

ggsave("paper/four_pndemics.png", dpi = 300, width = 297, height = 210, units = "mm")

Percent

plot_years <- c(seq(1890 - 1, 1890 + 2),
                seq(1918 - 1, 1918 + 2),
                seq(1957 - 1, 1957 + 2),
                seq(2020 - 1, 2020 + 2))

plot_data <- expected_deaths_monthly %>% 
  filter(Year %in% plot_years) %>% 
  mutate(Percent = ((Deaths - fit_flu_hc) / fit_flu_hc),
         Difference = ifelse(Deaths > fit_flu_hc, "More", "Less")) %>% 
  group_by(Country, Pandemic) %>% 
  mutate(time = row_number()) %>% 
  ungroup()

plot_data %>% 
  ggplot(aes(x = time)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_vline(xintercept = 13, colour = "grey80") +
  geom_vline(xintercept = 24, colour = "grey80") +
  geom_rect(xmin = as.Date(paste0(year, "-01-01")),
            xmax = as.Date(paste0(year, "-12-31")),
            ymin = -Inf, ymax = Inf,
            fill = "#ffffbf", alpha = 0.01, inherit.aes = FALSE) +
  geom_col(aes(y = Percent, fill = Difference), size = 0.5) +
  facet_grid(vars(Country), vars(Pandemic), scales = "free_x") + 
  scale_x_continuous(limits = c(min(plot_data$time)-1, max(plot_data$time)+1),
                     breaks = c(1, 13, 24, 36, 48), 
                     labels = c("-12m", "Start", "End", "+12m", "+24m")) +
  xlab("") + 
  ylab("% above predicted Deaths") +
  ggtitle("Excess Deaths as % of observed") +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(values = c("#67A9CF", "#EF8A62")) +  
  guides(fill = guide_legend(reverse = TRUE)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(), 
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(), 
        panel.grid.minor.y = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "grey40")) 

ggsave("paper/four_pndemics_perc.png", dpi = 300, width = 297, height = 210, units = "mm")

Table 1

pandemic_year <- c("1890","1918","1957","2020")

table_1 <- expected_deaths_monthly %>%
  select(Country,Year, Deaths, Population, fit_flu_hc) %>%
  filter(Year %in% pandemic_year)%>%
  group_by(Year, Country) %>%
  summarise(Expected_Mortality = round(sum(fit_flu_hc), 0),
            Deaths_num = sum(Deaths),
            Cum_excess_death = round(sum(Deaths)-sum(fit_flu_hc), 0),
            Percent_excess_death = round((Cum_excess_death/Expected_Mortality)*100, 1))%>%
  ungroup()

openxlsx::write.xlsx(table_1,"paper/table1.xlsx", row.names = FALSE, overwrite = TRUE)
Year Country Expected_Mortality Deaths_num Cum_excess_death Percent_excess_death
1890 Sweden 73077 81824 8747 12.0
1890 Switzerland 58397 61805 3408 5.8
1918 Spain 451779 695758 243979 54.0
1918 Sweden 78642 104591 25949 33.0
1918 Switzerland 50473 75034 24561 48.7
1957 Spain 282888 293502 10614 3.8
1957 Sweden 70064 73132 3068 4.4
1957 Switzerland 52280 51066 -1214 -2.3
2020 Spain 425484 479760 54276 12.8
2020 Sweden 88975 97870 8895 10.0
2020 Switzerland 67445 75570 8125 12.0

Appendix

## Crude
expected_deaths_monthly %>% 
  ggplot(aes(x = Date, y = Deaths)) +
  geom_col() +
  geom_rect(xmin = as.Date(paste0(1890, "-01-01")),
            xmax = as.Date(paste0(1890, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  geom_rect(xmin = as.Date(paste0(1918, "-01-01")),
            xmax = as.Date(paste0(1918, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  geom_rect(xmin = as.Date(paste0(1957, "-01-01")),
            xmax = as.Date(paste0(1957, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  geom_rect(xmin = as.Date(paste0(2020, "-01-01")),
            xmax = as.Date(paste0(2020, "-12-31")),
            ymin = 0, ymax = Inf,
            fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
  facet_grid(vars(Country), scales = "free_y") + 
  scale_y_continuous(labels = label_number(accuracy = 1L, suffix = "K", scale = 1e-4)) +
  theme_minimal() +
  xlab("") + ylab("Deaths") +
  ggtitle("Monthly number of Deaths")
## Per pop
expected_deaths_monthly %>% 
  ggplot(aes(x = Date, y = (Deaths/pop)*10000, colour = Country)) +
  geom_line() +
  theme_minimal() +
  xlab("") + ylab("Deaths / 10,000 population") +
  ggtitle("Monthly number of Deaths")

Deaths per population

Country panel

Switzerland

Impact of extremes

Using 1918 pandemic as an example of impact on estimates from 1919 - 1923.
Estimates with extreme year (red) and without (blue).
Converge after 5 year smoothing window stops having effect (ie. in 1924).

Figure

Table

Using 1918 pandemic as an example of impact on estimates from 1919.

Including 1918
Excluding 1918
Country Month Observed Predicted Difference Predicted Difference
Switzerland 1 5599 6602 1003 4905 -694
Switzerland 2 5324 6516 1192 5269 -55
Switzerland 3 5404 6467 1063 5293 -111
Switzerland 4 5451 6291 840 4965 -486
Switzerland 5 4976 5851 875 4454 -522
Switzerland 6 4166 5278 1112 3963 -203
Switzerland 7 3851 4878 1027 3610 -241
Switzerland 8 3674 4875 1201 3428 -246
Switzerland 9 3370 5295 1925 3416 46
Switzerland 10 4043 5948 1905 3575 -468
Switzerland 11 4465 6484 2019 3906 -559
Switzerland 12 4609 6669 2060 4385 -224
Sweden 1 10999 9006 -1993 7236 -3763
Sweden 2 8289 8407 118 7314 -975
Sweden 3 7781 8367 586 7176 -605
Sweden 4 8076 8524 448 6989 -1087
Sweden 5 7137 8236 1099 6725 -412
Sweden 6 5647 7396 1749 6276 629
Sweden 7 5590 6625 1035 5690 100
Sweden 8 5093 6547 1454 5205 112
Sweden 9 5083 7375 2292 5063 -20
Sweden 10 6120 8797 2677 5368 -752
Sweden 11 7289 9869 2580 6034 -1255
Sweden 12 7185 9803 2618 6773 -412
Spain 1 45826 53570 7744 41064 -4762
Spain 2 39277 51033 11756 41176 1899
Spain 3 44287 50970 6683 38800 -5487
Spain 4 38960 50583 11623 35987 -2973
Spain 5 37423 48628 11205 34907 -2516
Spain 6 36609 46998 10389 36063 -546
Spain 7 41756 48642 6886 38159 -3597
Spain 8 40557 54903 14346 39118 -1439
Spain 9 37199 63712 26513 38205 1006
Spain 10 37952 69199 31247 36880 -1072
Spain 11 37904 66942 29038 36988 -916
Spain 12 45002 59761 14759 38903 -6099

Impact of fit vs CI

Fit vs. LPI vs. UPI

Fit
LPI
UPI
Year Observed Predicted Difference Predicted Difference Predicted Difference
1890 61805 58397 -3408 54345 -7460 62582 777
1918 75034 50473 -24561 46630 -28404 54740 -20294
1957 51066 52280 1214 48119 -2947 57067 6001
2009 62476 61814 -662 58550 -3926 65094 2618
2020 75570 67445 -8125 63365 -12205 71772 -3798

2015 & 2020

Comparing BfS and our method for two years.

Fit
LPI
UPI
Year Method Observed Predicted Difference Predicted Difference Predicted Difference
2010 BfS 62038 61342 -696 55301 -6737 67387 5349
2010 Our estimates 62649 62476 -173 58927 -3722 66280 3631
2015 BfS 68065 65602 -2463 59369 -8696 71843 3778
2015 Our estimates 67606 65253 -2353 61814 -5792 68872 1266
2020 BfS 76391 68441 -7950 62117 -14274 74760 -1631
2020 Our estimates 75570 67445 -8125 63365 -12205 71772 -3798

Alt Viz

Stacked

Diffs 2